#pragma once
#include "../common.hpp"
#include "../Utility/Log.hpp"
#include "../Utility/IOObject.hpp"
#include "../Utility/Tools.hpp"
#include "Math.hpp"
#undef min
#undef max


namespace zzz{
template <unsigned int N, typename T>
class VectorBase
{
protected:
  T v[N];

public:
  // STL support
  // type definitions
  typedef T               value_type;
  typedef T*              iterator;
  typedef const T*        const_iterator;
  typedef T&              reference;
  typedef const T&        const_reference;
  typedef zsize           size_type;
  typedef size_type       difference_type;

  // iterator support
  iterator begin() { return v; }
  const_iterator begin() const { return v; }
  iterator end() { return v+N; }
  const_iterator end() const { return v+N; }

  // reverse iterator support
  typedef std::reverse_iterator<iterator> reverse_iterator;
  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;

  reverse_iterator rbegin() { return reverse_iterator(end()); }
  const_reverse_iterator rbegin() const { return const_reverse_iterator(end());}
  reverse_iterator rend() { return reverse_iterator(begin()); }
  const_reverse_iterator rend() const {return const_reverse_iterator(begin());}

  // at()
  reference at(size_type i) {
    ZRCHECK_LT(i, N);
    return v[i]; 
  }
  const_reference at(size_type i) const {
    ZRCHECK_LT(i, N);
  return v[i]; 
  }

  // front() and back()
  reference front() { return v[0]; }
  const_reference front() const {  return v[0];}
  reference back() {   return v[N-1]; }
  const_reference back() const {   return v[N-1]; }

  // size is constant
  static size_type size() { return N; }
  static bool empty() { return false; }
  static size_type max_size() { return N; }
  enum { static_size = N };

  // swap (note: linear complexity)
  void swap (VectorBase<N, T>& y) {  std::swap_ranges(begin(),end(), y.begin());}

  // check range (may be private because it is static)
  bool CheckIndex (size_type i) { return (i>=0 && i<size());}

  //subscript
  inline T& operator [](int a) {return at(a);}
  inline const T& operator [](int a) const {return at(a);}

  //copy
  inline const VectorBase<N, T> &operator=(const T &a){for (unsigned int i=0; i<N; ++i) v[i]=a; return *this;}
  inline const VectorBase<N, T> &operator=(const VectorBase<N, T>& a){IOObject<T>::CopyData(v, a.v, N); return *this;}
  template <typename T1>
  inline const VectorBase<N, T> &operator=(const VectorBase<N,T1>& a){for (unsigned int i=0; i<N; ++i) v[i]=(T)a[i]; return *this;}

  //raw data access
  inline T* Data(){return v;}
  inline const T* Data() const {return v;}

  //rawcopy
  inline void RawCopy(const VectorBase<N, T>& a){memcpy(v, a.v,sizeof(T)*N);}
  inline void RawCopy(const T *a){memcpy(v, a.v,sizeof(T)*N);}
  inline void Zero(zuchar x=0){memset(v, x, sizeof(T)*N);}
  inline static VectorBase<N, T> GetZero(int x=0){VectorBase<N, T> v; v.Zero(x); return v;}
  inline void Fill(const T &a){for (size_type i=0; i<N; ++i) v[i]=a;}
  inline static VectorBase<N, T> GetFill(const T &x){VectorBase<N, T> v; v.Fill(x); return v;}

  //positive and negative
  inline const VectorBase<N, T> operator+() const {return *this;}
  inline const VectorBase<N, T> operator-() const {VectorBase<N, T> ret; for (unsigned int i=0; i<N; ++i) ret.v[i]=-v[i]; return ret;}

  //addition
  inline const VectorBase<N, T> operator+(const VectorBase<N, T> &a) const {
    VectorBase<N, T> ret;
    for (unsigned int i=0; i<N; ++i)  ret.v[i]=v[i]+a.v[i];
    return ret;
  }
  inline const VectorBase<N, T> operator+(const T &a) const {
    VectorBase<N, T> ret;
    for (unsigned int i=0; i<N; ++i)  ret.v[i]=v[i]+a;
    return ret;
  }
  inline void operator+=(const VectorBase<N, T> &a) {
    for (unsigned int i=0; i<N; ++i)  v[i]+=a.v[i];
  }
  inline void operator+=(const T &a) {
    for (unsigned int i=0; i<N; ++i)  v[i]+=a;
  }

  //subtraction
  inline const VectorBase<N, T> operator-(const VectorBase<N, T> &a) const {
    VectorBase<N, T> ret;
    for (unsigned int i=0; i<N; ++i)  ret.v[i]=v[i]-a.v[i];
    return ret;
  }
  inline const VectorBase<N, T> operator-(const T &a) const {
    VectorBase<N, T> ret;
    for (unsigned int i=0; i<N; ++i)  ret.v[i]=v[i]-a;
    return ret;
  }
  inline void operator-=(const VectorBase<N, T> &a) {
    for (unsigned int i=0; i<N; ++i)  v[i]-=a.v[i];
  }
  inline void operator-=(const T &a) {
    for (unsigned int i=0; i<N; ++i)  v[i]-=a;
  }

  //multiplication
  inline const VectorBase<N, T> operator*(const VectorBase<N, T> &a) const {
    VectorBase<N, T> ret;
    for (unsigned int i=0; i<N; ++i)  ret.v[i]=v[i]*a.v[i];
    return ret;
  }
  inline const VectorBase<N, T> operator*(const T &a) const {
    VectorBase<N, T> ret;
    for (unsigned int i=0; i<N; ++i)  ret.v[i]=v[i]*a;
    return ret;
  }
  inline void operator*=(const VectorBase<N, T> &a) {
    for (unsigned int i=0; i<N; ++i)  v[i]*=a.v[i];
  }
  inline void operator*=(const T &a) {
    for (unsigned int i=0; i<N; ++i)  v[i]*=a;
  }

  //division
  inline const VectorBase<N, T> operator/(const VectorBase<N, T> &a) const {
    VectorBase<N, T> ret;
    for (unsigned int i=0; i<N; ++i)  ret.v[i]=v[i]/a.v[i];
    return ret;
  }
  inline const VectorBase<N, T> operator/(const T &a) const {
    VectorBase<N, T> ret;
    for (unsigned int i=0; i<N; ++i)  ret.v[i]=v[i]/a;
    return ret;
  }
  inline void operator/=(const VectorBase<N, T> &a) {
    for (unsigned int i=0; i<N; ++i)  v[i]/=a.v[i];
  }
  inline void operator/=(const T &a) {
    for (unsigned int i=0; i<N; ++i)  v[i]/=a;
  }

  //math
  inline T Dot(const VectorBase<N, T> &a) const {T ret=0;for (unsigned int i=0; i<N; ++i) ret+=v[i]*a.v[i]; return ret;}
  inline T Sum() const{T ret=0;for (unsigned int i=0; i<N; ++i)  ret+=v[i];return ret;}
  inline void Negative() {for (size_type i=0; i<N; ++i) v[i]=-v[i];}
  inline T Len() const {return Sqrt<T>((double)LenSqr());}
  inline T LenSqr() const {return FastDot(*this, *this);}

  inline T SafeNormalize() {T norm=Len(); if (norm!=0) *this/=norm; return norm;}
  inline VectorBase<N, T> SafeNormalized() const {T norm=Len(); if (norm!=0) return *this/norm; return *this;}
  inline T Normalize() {T norm=Len(); ZCHECK_NE(norm, 0)<<"Vector Normalize Failed!"; *this/=norm; return norm;}
  inline VectorBase<N, T> Normalized() const {T norm=Len(); ZCHECK_NE(norm, 0)<<"Vector Normalized Failed!"; return *this/norm;}

  inline T DistTo(const VectorBase<N, T> &a) const {return (T)sqrt((double)DistToSqr(a));}
  inline T DistToSqr(const VectorBase<N, T> &a) const {VectorBase<N, T> diff(*this-a); return diff.Dot(diff);}

  VectorBase<N, T> RandomPerpendicularUnitVec() const {
    VectorBase<N, T> ret((T)0);
    ret[Abs(v[0])>(T)(1-EPSILON)?1:0]=(T)1;
    ret-=*this*ret.Dot(*this)/LenSqr();
    ret.Normalize();
    return ret;
  }

  //Min and Max
  inline T Max() const {
    T maxv=v[0];
    for (unsigned int i=1; i<N; ++i)  if (maxv<v[i]) maxv=v[i];
    return maxv;
  }
  inline int MaxPos() const {
    T maxv=v[0];
    int pos=0;
    for (unsigned int i=1; i<N; ++i) if (maxv<v[i]) {maxv=v[i];pos=i;}
    return pos;
  }
  inline T Min() const {
    T minv=v[0];
    for (unsigned int i=1; i<N; ++i) if (minv>v[i]) minv=v[i];
    return minv;
  }
  inline int MinPos() const {
    T minv=v[0];
    int pos=0;
    for (unsigned int i=1; i<N; ++i) if (minv>v[i]) {minv=v[i];pos=i;}
    return pos;
  }
  inline T AbsMax() const {
    T maxv=Abs(v[0]);
    for (unsigned int i=1; i<N; ++i)  if (maxv<Abs(v[i])) maxv=Abs(v[i]);
    return maxv;
  }
  inline int AbsMaxPos() const {
    T maxv=Abs(v[0]);
    int pos=0;
    for (unsigned int i=1; i<N; ++i) if (maxv<Abs(v[i])) {maxv=Abs(v[i]);pos=i;}
    return pos;
  }
  inline T AbsMin() const {
    T minv=Abs(v[0]);
    for (unsigned int i=1; i<N; ++i) if (minv>Abs(v[i])) minv=Abs(v[i]);
    return minv;
  }
  inline int AbsMinPos() const {
    T minv=Abs(v[0]);
    int pos=0;
    for (unsigned int i=1; i<N; ++i) if (minv>Abs(v[i])) {minv=Abs(v[i]);pos=i;}
    return pos;
  }
  inline VectorBase<N, T> RotateToLeft() const {
    VectorBase<N, T> ret;
    for (size_type i = 0; i < N-1; ++i) ret[i] = v[i+1];
    ret[N-1] = v[0];
    return ret;
  }
  inline VectorBase<N, T> RotateToRight() const {
    VectorBase<N, T> ret;
    T tmp = v[N - 1];
    for (size_type i = N-1; i >= 1; --i) ret[i]=v[i-1];
    ret[0] = v[N - 1];
    return ret
  }
  inline VectorBase<N, T> Reverse() const {
    VectorBase<N, T> ret;
    for (size_type i = 0; i < N; ++i) ret[i] = v[N - 1 - i];
    return ret;
  }
  inline void KeepMin(const VectorBase<N, T> &y)
  {
    for (size_type i=0; i<N; ++i) v[i]=min(v[i], y[i]);
  }
  inline void KeepMax(const VectorBase<N, T> &y)
  {
    for (size_type i=0; i<N; ++i) v[i]=max(v[i], y[i]);
  }
  // comparisons
  bool operator== (const VectorBase<N, T>& y) const {return SafeEqual(begin(), end(), y.begin(), y.end());}
  bool operator< (const VectorBase<N, T>& y) const {return std::lexicographical_compare(begin(),end(), y.begin(), y.end());}
  bool operator!= (const VectorBase<N, T>& y) const {return !(*this==y);}
  bool operator> (const VectorBase<N, T>& y) const {return y<*this;}
  bool operator<= (const VectorBase<N, T>& y) const {return !(y<*this);}
  bool operator>= (const VectorBase<N, T>& y) const {return !(*this<y);}  
  
  //IO
  friend inline ostream& operator<<(ostream& os,const VectorBase<N, T> &me)
  {
    for (int i=0; i<N; ++i)
      os<<me.v[i]<<' ';
    return os;
  }
  friend inline istream& operator>>(istream& is,VectorBase<N, T> &me)
  {
    for (int i=0; i<N; ++i)
      is>>me.v[i];
    return is;
  }

  //constructor
  VectorBase(void){}
  VectorBase(const VectorBase<N, T>& a) {*this=a;}
  explicit VectorBase(const T &a){*this=a;}
  template<typename T1>
  explicit VectorBase(const T1 *p) {for (size_type i=0; i<N; ++i) v[i]=p[i];}
  template<typename T1>
  explicit VectorBase(const VectorBase<N,T1>& a) {*this=a;}
  template<typename T1>
  explicit VectorBase(const VectorBase<N-1,T1>& a, const T1 &b) {for (size_type i=0; i<N-1; ++i) v[i]=a[i]; v[N-1]=b;}
  template<typename T1>
  explicit VectorBase(const VectorBase<N+1,T1>& a)  {for (size_type i=0; i<N; ++i) v[i]=a[i];}
};

//scale at front
template <zuint N, typename T>
inline const VectorBase<N, T> operator+(const T &v,const VectorBase<N, T> &vec)
{
  return vec+v;
}
template <zuint N,typename T>
inline const VectorBase<N, T> operator-(const T &v,const VectorBase<N, T> &vec)
{
  return -vec+v;
}
template <zuint N,typename T>
inline const VectorBase<N, T> operator*(const T &v,const VectorBase<N, T> &vec)
{
  return vec*v;
}
template <zuint N,typename T>
inline const VectorBase<N, T> operator/(const T &v,const VectorBase<N, T> &vec)
{
  VectorBase<N, T> res;
  for (size_type i=0; i<N; ++i) res[i]=v/vec[i];
  return res;
}



//write all these stuff just to avoid direct memory copy when construct and copy
template <unsigned int N,typename T>
class Vector : public VectorBase<N, T>
{
protected:
  using VectorBase<N, T>::v;
public:
  //constructor
  Vector(void){}
  Vector(const Vector<N, T>& a):VectorBase<N, T>(a) {}
  Vector(const VectorBase<N, T>& a):VectorBase<N, T>(a) {}
  explicit Vector(const T &a):VectorBase<N, T>(a){}
  explicit Vector(const VectorBase<N-1, T> &a, const T &b):VectorBase<N, T>(a, b){}
  explicit Vector(const VectorBase<N+1, T> &a):VectorBase<N, T>(a){}
  template<typename T1> explicit Vector(const T1 *p):VectorBase<N, T>(p) {}
  template<typename T1> explicit Vector(const Vector<N,T1>& a):VectorBase<N, T>(a) {}
  template<typename T1> explicit Vector(const VectorBase<N,T1>& a):VectorBase<N, T>(a) {}

  //assign
  inline const Vector<N, T> &operator=(const Vector<N, T>& a){VectorBase<N, T>::operator =(a); return *this;}
  using VectorBase<N, T>::operator=;

};

//meta-programming...
// primary template
template <unsigned int N, typename T>
struct DotProductHelper 
{
  static T Result (const T* a, const T* b) 
  {
    return *a * *b  +  DotProductHelper<N-1, T>::Result(a+1, b+1);
  }
};

// partial specialization as end criteria
template <typename T>
struct DotProductHelper<1, T> 
{
  static T Result (const T* a, const T* b) {return *a * *b;}
};

// convenience function
template <unsigned int N, typename T>
inline T FastDot(const VectorBase<N, T> &a, const VectorBase<N, T> &b)
{
  return DotProductHelper<N, T>::Result(a.Data(), b.Data());
}

// IOObject

template<unsigned int N, typename T>
class IOObject<VectorBase<N, T> >
{
public:
  static void WriteFileB(FILE *fp, const VectorBase<N, T>* src, zsize size) {
    IOObject<T>::WriteFileB(fp, src->Data(), size*N);
  }
  static void ReadFileB(FILE *fp, VectorBase<N, T>* dst, zsize size) {
    IOObject<T>::ReadFileB(fp, dst->Data(), size*N);
  }
  static void WriteFileB(FILE *fp, const VectorBase<N, T>& src) {
    IOObject<T>::WriteFileB(fp, src.Data(), N);
  }
  static void ReadFileB(FILE *fp, VectorBase<N, T>& dst) {
    IOObject<T>::ReadFileB(fp, dst.Data(), N);
  }

  static void WriteFileR(RecordFile &fp, const zint32 label, const VectorBase<N, T>& src) {
    IOObject<T>::WriteFileR(fp, label, src.Data(), N);
  }
  static void ReadFileR(RecordFile &fp, const zint32 label, VectorBase<N, T>& dst) {
    IOObject<T>::ReadFileR(fp, label, dst.Data(), N);
  }
  static void WriteFileR(RecordFile &fp, const VectorBase<N, T>& src) {
    IOObject<T>::WriteFileR(fp, src.Data(), N);
  }
  static void ReadFileR(RecordFile &fp, VectorBase<N, T>& dst) {
    IOObject<T>::ReadFileR(fp, dst.Data(), N);
  }

  static void WriteFileR(RecordFile &fp, const zint32 label, const VectorBase<N, T>* src, zsize size) {
    IOObject<T>::WriteFileR(fp, label, src->Data(), size*N);
  }
  static void ReadFileR(RecordFile &fp, const zint32 label, VectorBase<N, T>* dst, zsize size) {
    IOObject<T>::ReadFileR(fp, label, dst->Data(), size*N);
  }
  static void WriteFileR(RecordFile &fp, const VectorBase<N, T>* src, zsize size) {
    IOObject<T>::WriteFileR(fp, src->Data(), size*N);
  }
  static void ReadFileR(RecordFile &fp, VectorBase<N, T>* dst, zsize size) {
    IOObject<T>::ReadFileR(fp, dst->Data(), size*N);
  }
  static void CopyData(VectorBase<N, T>* dst, const VectorBase<N, T>* src, zsize size) {
    IOObject<T>::CopyData(static_cast<T*>(dst->Data()), static_cast<const T*>(src->Data()), size*N);
  }
};

template<unsigned int N, typename T>
class IOObject<Vector<N, T> >
{
public:
  static void WriteFileB(FILE *fp, const Vector<N, T>* src, zsize size)
  {
    IOObject<VectorBase<N, T> >::WriteFileB(fp, src, size);
  }
  static void ReadFileB(FILE *fp, Vector<N, T>* dst, zsize size)
  {
    IOObject<VectorBase<N, T> >::ReadFileB(fp, dst, size);
  }
  static void WriteFileB(FILE *fp, const Vector<N, T>& src)
  {
    IOObject<VectorBase<N, T> >::WriteFileB(fp, src);
  }
  static void ReadFileB(FILE *fp, Vector<N, T>& dst)
  {
    IOObject<VectorBase<N, T> >::ReadFileB(fp, dst);
  }

  static void WriteFileR(RecordFile &fp, const zint32 label, const Vector<N, T>& src) {
    IOObject<VectorBase<N, T> >::WriteFileR(fp, label, src);
  }
  static void ReadFileR(RecordFile &fp, const zint32 label, Vector<N, T>& dst) {
    IOObject<VectorBase<N, T> >::ReadFileR(fp, label, dst);
  }
  static void WriteFileR(RecordFile &fp, const Vector<N, T>& src) {
    IOObject<VectorBase<N, T> >::WriteFileR(fp, src);
  }
  static void ReadFileR(RecordFile &fp, Vector<N, T>& dst) {
    IOObject<VectorBase<N, T> >::ReadFileR(fp, dst);
  }

  static void WriteFileR(RecordFile &fp, const zint32 label, const Vector<N, T>* src, zsize size) {
    IOObject<VectorBase<N, T> >::WriteFileR(fp, label, src, size);
  }
  static void ReadFileR(RecordFile &fp, const zint32 label, Vector<N, T>* dst, zsize size) {
    IOObject<VectorBase<N, T> >::ReadFileR(fp, label, dst, size);
  }
  static void WriteFileR(RecordFile &fp, const Vector<N, T>* src, zsize size) {
    IOObject<VectorBase<N, T> >::WriteFileR(fp, src, size);
  }
  static void ReadFileR(RecordFile &fp, Vector<N, T>* dst, zsize size) {
    IOObject<VectorBase<N, T> >::ReadFileR(fp, dst, size);
  }
  static void CopyData(Vector<N, T>* dst, const Vector<N, T>* src, zsize size)
  {
    IOObject<VectorBase<N, T> >::CopyData(dst, src, size);
  }
};

template <zuint N, typename T> inline Vector<N, T> Abs(const VectorBase<N, T> &x) {
  Vector<N, T> ret;
  for (size_type i=0; i<N; ++i) ret[i]=Abs(x[i]);
  return ret;
}

}
